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Using a linear analysis, we study the stability of giant-vortex states in very thin 
disks. The vortex expulsion and penetration fields are obtained for finite thickness 
disks from a numerical solution of the non-linear Ginzburg-Landau (GL) equations. 
Using an extension of the London approximation, in which the phase distribution of 
the order parameter is prescribed and the superconducting electron density is found 
numerically, we consider the free energy behavior for transitions between different 
i Q | superconducting states. 
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Recently, much attention is devoted to investigations of superconducting phe- 
nomena in mesoscopic disks whose radius R and thickness d are comparable 
to the coherence £ and penetration A lengths. Precise measurements of the 
magnetization of single superconducting Al disks using Hall-magnetometry 
allowed to follow the evolution of few vortex states [1,2]. A number of re- 
markable observations (the influence of disk size on the type and order of the 
phase transition between the normal and the superconducting state, promi- 
nent hysteretic behavior in defect-free samples, the paramagnetic Meissner 
effect in the field-cooling regime, etc.) initiated several theoretical works, in 
which superconductivity of mesoscopic samples is considered within the full 
GL approach. The superconduncting state may present either a giant vortex 
state with an axially symmetric distribution of the superconducting density 
or a multi-vortex Abrikosov-like state, which is usually realized in type-// su- 
perconductors. Since the effective penetration length A* = X 2 /d [3] may by far 
exceed that of a bulk superconductor, the multi-vortex states can also exist 
in thin disks made from type-/ superconducting material. 

The giant vortex states for finite thickness disks were considered in Refs. [4,5]. 
With increasing disk radius a second-order reversible phase transition from 
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the normal to the superconducting state observed for small disk radii is re- 
placed by first-order transitions with jumps in the magnetization [4,5]. The 
simulated magnetizations are in good quantitative agreement with those from 
experiment [4,6,7]. The multi-vortex states have been investigated using the 
London approach [8], the GL approach [6,7,9], and the lowest Landau level 
approximation [10]. The free energy barriers separating the superconducting 
states with different vorticity were obtained using an expansion of the order 
parameter over the eigenfunctions of the GL kinetic energy operator [11]. The 
lowest Landau level approximation was used to study these barriers in Ref.[12]. 
An equilibrium phase diagram, showing which state is energetically preferable 
as function of the magnetic field and disk thickness, was found in Ref.[9]. In 
the present work, we study a non- equilibrium phase diagram showing the sta- 
bility region of the giant vortex states as well as the order of the transition 
between different superconducting states. 

We consider a superconducting disk surrounded by an insulator medium and 
placed in a uniform magnetic field H, which is perpendicular to the disk 
surface. Measuring the distance in units of the coherence length £, the vector 



potential A in ch/2e^, and the order parameter \1/ in J— a/ ' (3 with a, (3 being 
the GL coefficients [3], we write the system of GL equations for a thin (d <C 
£, A) disk in the following form [9]: 
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where the indices 2D, 3D refer to two-dimensional (in the disk plane (x,y)) 
and three-dimensional operators; k, = A/£ is the GL parameter; j 2 o is the 
density of superconducting current. The boundary conditions to Eqs. (1,2) 
correspond to zero superconducting current at the disk boundaries and an 
uniform external magnetic field far from the disk A^^^ — |[r X Ho], respec- 
tively. Our numerical approach to solve Eqs. (1,2) is described in [9]. 

We restrict our considerations of the stability of the giant vortex states to 
the case of small disks A* ^> R, when the expulsion of the magnetic field 
can be neglected. Then the order parameter of the giant vortex state can be 
presented as * = ifj(p)e lL ^, where L is the angular quantum momentum; p, <fi 
are the cylindrical coordinates. The radial wavefunction ip obeys the following 
equation 

I d dip .L 1 TT , 2 , ,, 

—ir pit + ~ + o h °p ^ = ^ - ^ > 4 

pop op p 2 



with the boundary condition (dip/dp) p= R = 0. In order to investigate whether 
the giant vortex state is stable with respect to small perturbations from other 
angular momentum states we use the time-dependent (TD) first GL equation 
[13]. Representing the order parameter as a mixture of three angular harmon- 
ics * = i){p)e iL ^ + Si(p)e iLl<f,+flt + S 2 {p)e iL ^ + ^, we linearize the first TDGL 
equation with respect to small perturbations Si, 62 and obtain the following 
set of equations 
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which describes our eigen-value problem. Here, 

Gl = -~p^- + (- + \h qP )<2 

pop op p I 

is the Schrodinger operator in an uniform magnetic field with the boundary 
condition {dSi^ldp)p=R = 0. The sign of the eigenvalue p determines whether 
the considered state is stable p < or unstable p > 0. Note, that the use 
of three harmonics and the relation Li + L 2 = 2L between them is dictated 
by the non-linear term in the RHS of Eq.(l), which corresponds to the third 
power of the order parameter. Generally, a linear analysis does not allow to 
predict which state will be realized from the evolution of instability. However, 
according to our previous results from solving the non-linear GL equations 
[9] we expect that the perturbation with L x = 0, L 2 = 1L leaves the vor- 
ticity unchanged and leads to the appearance of a state corresponding to L 
vortices arranged in a ring. When the angular momentum of the perturba- 
tion is not a multiple of that of the initial state, we expect the appearance 
of a supercondicting state with a different vorticity. To find the perturbation 
spectrum we use a finite-difference representation of the Schrodinger operator 
G [5] and reduce Eq.(5) to a matrix, which is numerically diagonalized with 
the Householder technique. The state becomes unstable if the maximum value 
of p found for different pairs Li,L2 = 1L — Li changes its sign and becomes 
positive. 

The results from our linear analysis are shown in Fig. 1 for L = 0, ...,9. The 
Meissner state becomes unstable relative to the entrance of extra vortices with 
increasing magnetic field above the penetration field H p . For small disk radii, 
it is the perturbation with Li — L — 1,L 2 — L + 1 which makes the L state 
unstable. With further increasing the disk radius, our linear analysis shows 
that a couple of vortices will enter the system at once. The points, at which 
the number of penetrating vortex changes, are denoted by solid squares in 
Fig. 1. A similar behavior of vortex penetration is also observed for the single 
vortex state. It should be noted, that the eigenvalues of the perturbations 



with different angular momenta are very close to each other at large disk radii. 
Therefore, a non-linear consideration is required to answer how many vortices 
can simultaneously enter into the disk. Unfortunately, such an analysis is a 
very difficult task since the results turn out to be very sensitive to the shape 
of the perturbation. Our numerical solution of the non-linear GL equations 
shows that occasionally several vortices enter indeed into the system. Although 
we are able to predict the entrance field with high accuracy, the number of 
penetrating vortices is rather uncertain. This indicates the at the entrance 
field the barriers separating the different L-states become very small. 

The Meissner state and the single vortex state can exist in arbitrary large 
disks. The expulsion field H e , below which the vortex state becomes unstable 
relative to the transition to the Meissner state, decreases with disk radius, 
in such a way that the external magnetic flux piercing the disk tends to the 
flux quantum $ — hc/2e. The other giant vortex states remain stable only 
inside some region of the phase diagram. Their stability is restricted by either 
the transition to the multi-vortex state with the same vorticity, which occurs 
as a rule with decreasing magnetic field, or the transition to another giant 
vortex state. The left (right) boundary of the stability region corresponds to 
expulsion (penetration) of a single vortex with decreasing (increasing) disk 
radius. Note, that the state with L = 2 in principle allows for a reentrant 
behavior - the transition to the multi-vortex state with both decreasing and 
increasing magnetic field. Unfortunately, this prediction is not confirmed by 
our non-linear analysis. The boundary of transition between the giant and 
multi-vortex state increases with angular momentum and, as expected, tends 
to saturate at the second critical field H c i for large disk radii. The stability 
regions of different states overlap strongly and many different superconducting 
states can exist with the same external conditions. 

To study the influence of the disk thickness on the penetration and expulsion 
field we solve numerically Eqs. (1,2). Starting from the Meissner (vortex) state 
we increase (decrease) slowly the applied magnetic field. Using the supercon- 
ducting state obtained at the previous step as initial condition, we find a new 
steady-state solution of Eqs. (1,2) and check whether the vorticity changes. 
The penetration and expulsion fields found for zero disk thickness coincide 
with those from the linear analysis above. It turns out that the vortex state 
exibits a weak paramagnetic Meissner effect just before vortex expulsion. This 
is the main reason why the expulsion field (thin dotted curves in Fig. 2) de- 
creases with disk thickness. The Meissner state shows, as expected, a strong 
diamagnetic response. The decrease of the total magnetic flux piercing through 
the disk results in an expansion of the stability region of the Meissner state 
with increasing d/X 2 . Recall, that the field expulsion from thin disks is only 
determined by the ratio d/X 2 , which accounts for both the disk thickness and 
the penetration length of bulk material. The first critical field H c i, above which 
the vortex state becomes energetically more favorable, also increases with the 



disk thickness. All the critical fields H cl = 2H c2 t,/R, H e = 2H c2 t, 2 /R 2 , and 
H p = 2H C 2^ln(R/£ s )/R can be found within the London approximation in the 
limit of zero disk thickness [8,11]. Note, that the penetration field from the 
London approximation does not agree with that from the GL theory because 
the saddle point state for vortex entrance does not correspond to any vortex 
state nearby H p [11]. However, the London approach predicts rather accurately 
the expulsion and the first critical field for R ^> £ disk radii. 

The superconducting density is assumed to be a constant in the conventional 
London approximation. This assumption breaks down at small disk radii or 
small inter-vortex distances which is the case of mesoscopic samples. An ex- 
tension of the London approach accounting for variation of the amplitude ip of 
the order parameter ^(p) = il)(p)e %s ^ has been proposed in [11]. The phase 
distribution S, which depends on the vortex positions pi, is assumed to be cre- 
ated by the vortices and their mirror images with coordinates ft and piR 2 /p 2 , 
respectively [8,11]. Then the amplitude of the order parameter obeys the first 
GL equation 

-A-0 + (VS - A) 2 i) = i) - ^ 3 , (6) 

with the boundary condition (dip/dp) p= R = 0. Solving numerically the last 
equation we find the free energy F(pi) = — f ip dr/ndR 2 , where integration 
is performed over the disk volume. Doing so we reduce the GL free energy 
functional defined in the functional space to the functional of the vortex co- 
ordinates. Then the vortices can be considered as classical particles, whose 
motion is governed by this free energy. Such an approach, which we call the 
modified London approach below, is shown to lead to accurate results in the 
case when vortices do not cross the sample boundary [11]. 

According to numerical simulations [9], the transition from the giant vortex 
state to the multi-vortex state with the same vorticity is followed by a weak 
jump in the slope of the magnetization indicating a second-order transition. 
The modified London approach allows to study the free energy as function of 
the vortex coordinates, which is shown in Fig. 3 for the two-vortex state. At 
large magnetic fields, the minimum of the free energy is achieved when both 
vortices are located in the disk center. With decreasing magnetic field, the 
curvature of the potential curve decreases and tends to zero at some critical 
point H/H C 2 ~ 0.58. Thereafter, the multi-vortex state with separated vor- 
tices becomes more energetically favorable. No two minima of the free energy 
exist simultaneously. This agrees with results from our numerical solution of 
the GL equations [9] and proves that the transition is of second-order. 

At small disk radii, all the multi-vortexs state correspond to vortices arranged 
in a ring. With increasing disk radius, the configuration (1 : L — l) with a single 
vortex inside the ring also becomes stable at large vorticity and coexists with 



the ring structure (0 : L) [14] . There is experimental evidence of such a coex- 
istence [16]. Solving numerically Eqs. (1,2) we find the magnetization of both 
state (0 : 7), (1 : 6) (see Fig. 4). For R = 5£, the state with (0 : 7) turns out to 
be energetically preferable. Note, that the magnetization curves do not merge 
with each other suggesting a first-order transition between these multi-vortex 
states. To find the free energy barrier separating the states (1 : 6), (0 : 7) we 
apply the modified London approach. Starting from the configuration (1:6) 
we move slowly the central vortex to the disk boundary. The coordinates of 
the other vortices are found by minimizing the free energy. At small magnetic 
fields, there exists two minima of the free energy separated by a barrier (see 
Fig. 5). At some critical field H/H C 2 ~ 0.82 this barrier disappears and the 
state (1 : 6) becomes unstable. The transition between states (1 : 6) and (0 : 7) 
is followed by a jump in both the magnetization and the free energy which 
also proves the first-order of the transition between the multi-vortex states. 
Note, that the barrier separating two multi-vortex states is much smaller than 
that for vortex expulsion and penetration. 
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Fig. 1. The stability region of the different giant vortex states for zero disk thickness. 
The upper and lower solid thick curves correspond to the nucleation (H c ^) and 
expulsion (H e ) field, respectively. The thick-dashed curves denote the penetration 
fields (H p ) for the Meissner state and the single vortex state. The thin solid and 
dotted curves correspond to the transitions between the giant vortex states with 
different angular momenta and the transition from the giant vortex state to the 
multi-vortex state with the same vorticity (see inset) , respectively. 



Fig. 2. The first critical (H c i), penetration (H p ), and expulsion (H p ) field as function 
of the disk radius for different disk thicknesses. Thick curves denote results found 
with a conventional London approach. 



Fig. 3. The free energy of the two-vortex state as a function of the inter-vortex 
distance (a) for different magetic fields. 



Fig. 4. The unitless magnetization M = J(H — Ho)df> / Air 2 H&R? of the multivortex 
state L = 7 for two different vortex configurations and two thicknesses of the disk. 
The calculated magnetization includes the detector size effect [7] (square detector 
with width ~ 3.1//m). 



Fig. 5. The free energy as a function of the position of one of the vortices, which 
is shifted from the disk center to the disk boundary. The inset shows the vortex 
configurations in the two stable states. 
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